### This script will create graphs of climate data from multiple sites for multiple months
### Each site month will be produced as an independent graph to aid the visual checking of errors
### Script produced by Collin Storlie on 3rd November 2010 and updated on 15th February 2010

# First define the working directories

in.dir = "/home1/99/jc152199/underlog/"
setwd(in.dir)
out.dir = "/home1/99/jc152199/underlog/plots/experiment/"

# Import raw data, check the values, and examine the class of each column of the data frame

air.data = read.csv('/home1/99/jc152199/MicroclimateStatisticalDownscale/ToAnalyse/MicroMacroMinMaxASCII.csv', header=T)
raw.data = read.csv(paste(in.dir,'underlograwdata.csv',sep=''), header=T)
head(raw.data)
str(raw.data)

#Subset data to a few sites only for ease of processing

#raw.sub = subset(raw.data, site=='WU9A2')
#raw.sub2 = subset(raw.data, site=='WU11A2')
#raw.data = rbind(raw.sub, raw.sub2)

#Create columns for year, month, and day, populate them all with NA's

air.data$year = air.data$month = air.data$day = raw.data$datetime = NA

#Populate these new date columns by formatting the already existing complete datetime field first into a class date, then into character (with only the relevant portion of date)
#Then format these characters into numerals for easy sorting during the plotting process

air.data$year = as.numeric(format(as.Date(air.data$date,'%Y-%m-%d'),"%Y"))
air.data$month = as.numeric(format(as.Date(air.data$date,'%Y-%m-%d'),"%m"))
air.data$day = as.numeric(format(as.Date(air.data$date, '%Y-%m-%d'),"%d"))
air.data$airrange = air.data$micro_max - air.data$micro_min

raw.data$datetime = paste(raw.data$date, raw.data$time, sep=' ')

raw.data$datetimeformat = as.POSIXlt(raw.data$datetime, format= "%Y-%m-%d %H:%M")

#Aggregate readings based on site, year, month, day AND mission using the summary function to produce temperature range for under log data
# First define a function which produces the temperature range (i.e. max - min)

range.fun = function(x){return(max(x,na.rm=T)-min(x,na.rm=T))}

tdata.max = aggregate(raw.data$logtemp,by=list(site=raw.data$site, year=raw.data$year, month=raw.data$month, day=raw.data$day, mission=raw.data$mission), FUN = max)
tdata.min = aggregate(raw.data$logtemp,by=list(site=raw.data$site, year=raw.data$year, month=raw.data$month, day=raw.data$day, mission=raw.data$mission), FUN = min)
tdata.range = aggregate(raw.data$logtemp,by=list(site=raw.data$site, year=raw.data$year, month=raw.data$month, day=raw.data$day, mission=raw.data$mission), FUN = range.fun)

#Rename columns appropriately and then merge datasets

names(tdata.max)[6]='underlogmax'
names(tdata.range)[6]='ulrange'
names(tdata.min)[6]='underlogmin'

tdata = merge(tdata.min,tdata.max, by=c('site','year','month','day','mission'))
tdata = merge(tdata,tdata.range,by=c('site','year','month','day','mission'))

all.data = tdata

sites = unique(raw.data$site)

#Begin a loop that will plot climate data in chronological order for each site-month, combining all site-months for all years into a single .pdf document

for (xsite in sites) 

{

  sub.tdata = subset(all.data, site==xsite)
  years = unique(sub.tdata$year)
  missions = sort(unique(sub.tdata$mission), decreasing=F) 
  cols = seq(1:NROW(missions))
  mis.cols = cbind(missions,cols)
  mis.cols = as.data.frame(mis.cols, missions=mis.cols[,1], cols=mis.cols[,2])
  
  if (nrow(sub.tdata)>0) 
  
  {

  pdf(file=paste(out.dir,xsite,".pdf", sep=""))
  
  cat(xsite,'\n')

    for (xyear in years) 
    
    {
  
      sub2.tdata = subset(sub.tdata, year==xyear)
      
      month.list = sort(unique(sub2.tdata$month), decreasing=F)
      
      if (nrow(sub2.tdata)>0)
      
      cat(xyear,'\n')    
      
      {

        for (xmonth in month.list) 
        
        {
    
          sub3.tdata = subset(sub2.tdata, month==xmonth)

          if  (nrow(sub3.tdata)>1)                                                                                  
          
          {
          
          sm2 = data.frame(missions=NA,cols=NA)
          
            for (i in sort(unique(sub3.tdata$mission), decreasing=F))
          
            {
          
            sm2 = rbind(sm2,mis.cols[which(mis.cols$missions==i),])
                             
            }
          
          sm2 = sm2[2:nrow(sm2),]       
          plot(sub3.tdata$day, sub3.tdata$underlogmin, type='p', col=sm2[,2], pch = sm2[,2], ylim=c(min(sub3.tdata$underlogmin),max(sub3.tdata$underlogmin)+3), xlab = "Day", ylab = "Under Log Min Temp", main=paste(xsite," ",xyear,"-",xmonth,sep=""))
          legend(1,max(sub3.tdata$underlogmin)+2,c(sm2[,1]),pch=c(sm2[,2]),col=c(sm2[,2]))
          
          }
    
        }
    
      } 

    }

  }
                        
  dev.off()

}

# Bonus commands to put in the loop

#text(sub3.tdata$airrange, sub3.tdata$ulrange, labels=sub3.tdata$mission, cex=.5, pos=2)

#plot(sub3.tdata$day, sub3.tdata$underlogmax, type='n', col='black', xlab = "Day", ylab = "Under Log Max Temp", main=paste(xsite," ",xyear,"-",xmonth,sep=""), sub=paste('r2=',round(as.numeric(summary(lm1)[8]),3),' slope=',round(as.numeric(lm1$coefficients[2]),3),sep=''))

# A bonus loop to plot slope between airrange and ulrange for a moving window of 1 week

roundmean.fun = function(x){return(round(mean(x),0))}

#Begin a loop that will plot the relationship between log temp range and air temp range against time for each site individually

for (xsite in sites) 

{

  #sub.tdata = subset(tdata, site==xsite)
  #sub.airdata = subset(air.data, site==xsite)
  sub.tdata = subset(all.data, site==xsite)
  years = unique(sub.tdata$year)
  
  if (nrow(sub.tdata)>0) 
  
  {

  pdf(file=paste(out.dir,xsite,".pdf", sep=""))
  
  t.plotdata = data.frame(slope=NA, julavg=NA)
  
  cat(xsite,'\n')

    for (ii in seq(340,1126,7))
    
    {
    
    sub2.tdata = sub.tdata[which(sub.tdata$Julian>=ii & sub.tdata$Julian<=ii+27),]   # This command subsets the data based on Julian dates between ii and ii+27
    
      if (nrow(sub2.tdata)>1)
      
      {
    
      lm1 = lm(ulrange ~ airrange, data=sub2.tdata) # Find slope of relationship
    
      t.julavg = aggregate(sub2.tdata$Julian, by=list(site=sub2.tdata$site), FUN=roundmean.fun) # Calculate average julian date (proxy for time in weeks)
   
      slope = round(as.numeric(lm1$coefficients[2]),3)
    
      julavg = t.julavg[1,2]
    
      t.plotdata = rbind(t.plotdata,cbind(slope, julavg))
      
      }
    
    }
     
   }
   
   if (nrow(t.plotdata)>1)
   
   {
   
   plot(t.plotdata$julavg, t.plotdata$slope, type='b', col='black', xlab = "Juldate", ylab = "Slope", main=paste(xsite,sep=""))
   dev.off()
   
   }
   
  } 
